

#### Cross-national df construction #### 
#### gpfc 2020
#### 8-1-2020

#### Packages #### 

packs <- c('dplyr', 'lfe', 'ggplot2', 'ggthemes', 'cowplot')
lapply(packs, require, character.only = T)
setwd('~/Dropbox/public_health_repression/replication-for-IS/')
theme_set(theme_tufte())
df <- readRDS('data/cross_national_df.rds')


#### Descriptive Figures #### 
#### Examples of Countries #### 
setwd('~/Dropbox/public_health_repression/replication-for-IS/figures')


p1 <- ggplot(subset(df,cname=="South Africa"), aes(event_time, repress)) + 
  geom_line() +
  geom_vline(xintercept = 0, col='red') + 
  ggtitle("South Africa") +
  xlab("Months until COVID-19 Lockdown") +
  ylab("Repression Events") +
  theme_tufte()

p2 <- ggplot(subset(df,cname=="Nigeria"), aes(event_time, repress)) + 
  geom_line() +
  geom_vline(xintercept = 0, col='red') + 
  ggtitle("Nigeria") +
  xlab("Months until COVID-19 Lockdown") +
  ylab("Repression Events") +
  theme_tufte()

p3 <- ggplot(subset(df,cname=="Kenya"), aes(event_time, repress)) + 
  geom_line() +
  geom_vline(xintercept = 0, col='red') + 
  ggtitle("Kenya") +
  xlab("Months until COVID-19 Lockdown") +
  ylab("Repression Events") +
  theme_tufte()

p4 <- ggplot(subset(df,cname=="Uganda"), aes(event_time, repress)) + 
  geom_line() +
  geom_vline(xintercept = 0, col='red') + 
  ggtitle("Uganda") +
  xlab("Months until COVID-19 Lockdown") +
  ylab("Repression Events") +
  theme_tufte()

p5 <- cowplot::plot_grid(p1, p2, p3, p4, labels = "AUTO")
p5
ggsave("cn_examples.png", p5, width = 13)

plot_df <- df %>% 
  group_by(event_time) %>% 
  summarise(repress = sum(repress))

p <- ggplot(plot_df, aes(event_time, repress)) + 
  geom_line() +
  geom_vline(xintercept = 0, col='red') + 
  ggtitle("All African Countries in Sample") +
  xlab("Months Until COVID-19 Lockdowns") +
  ylab("Repression Events") +
  theme_tufte()

ggsave("cn_main.png", p, width = 13)


#### Main Results Figure #### 


m1 <- felm(repress ~ shutdown*event_time | factor(COWcode) | 0 | COWcode, 
           data = df)
m2 <- felm(repress ~ shutdown*event_time + I(case_lag/1000) + death_rate_lag + repress_lag | factor(COWcode) | 0 | COWcode, 
           data = df)
m3 <- felm(repress ~ shutdown*event_time + shutdown*I(event_time^2)  + I(case_lag/1000) + death_rate_lag + repress_lag| factor(COWcode) | 0 | COWcode, 
           data = df)
m4 <- felm(repress_bin ~ shutdown*event_time | factor(COWcode) | 0 | COWcode, 
           data = df)
m5 <- felm(repress_bin ~ shutdown*event_time + I(case_lag/1000) + death_rate_lag + repress_lag | factor(COWcode) | 0 | COWcode, 
           data = df)
m6 <- felm(repress_bin ~ shutdown*event_time + shutdown*I(event_time^2)  + I(case_lag/1000) + death_rate_lag + repress_lag | factor(COWcode) | 0 | COWcode, 
           data = df)

# Put model estimates into temporary data.frames:
model1Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m1$beta[1],
                          SE = m1$cse[1],
                          Models = "Baseline")
model2Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m2$beta[1],
                          SE = m2$cse[1],
                          Models = "Covariates")
model3Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m3$beta[1],
                          SE = m3$cse[1],
                          Models = "Quadradic")
# Combine these data.frames
allModelFrame <- data.frame(rbind(model1Frame, model2Frame, model3Frame))  # etc.

# Specify the width of your confidence intervals
interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

# Plot
zp1 <- ggplot(allModelFrame, aes(shape = Models))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
zp1 <- zp1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2), fill = "WHITE")
zp1 <- zp1 + theme_tufte()
zp1 <- zp1 + ggtitle("Repression (Levels)")  + xlab("") + ylab("Shutdown Coefficient") + 
  scale_fill_discrete(name = 'Model') + theme(legend.position = "bottom")
zp1



# Put model estimates into temporary data.frames:
model1Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m4$beta[1],
                          SE = m4$cse[1],
                          Models = "Baseline")
model2Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m5$beta[1],
                          SE = m5$cse[1],
                          Models = "Covariates")
model3Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m6$beta[1],
                          SE = m6$cse[1],
                          Models = "Quadradic")
# Combine these data.frames
allModelFrame <- data.frame(rbind(model1Frame, model2Frame, model3Frame))  # etc.

# Specify the width of your confidence intervals
interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

# Plot
zp2 <- ggplot(allModelFrame, aes(shape = Models))
zp2 <- zp2 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
zp2 <- zp2 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp2 <- zp2 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2), fill = "WHITE")
zp2 <- zp2 + theme_tufte()
zp2 <- zp2 + ggtitle("Repression (Binary)") + xlab("") + ylab("Shutdown Coefficient") + 
  scale_fill_discrete(name = 'Model') + theme(legend.position = "bottom")

plot_1 <- cowplot::plot_grid(zp1, zp2, labels = 'AUTO')

ggsave('base_plot.png', plot_1, width = 13)




#### Protest Results Figure #### 

m1 <- felm(protest ~ shutdown*event_time | factor(COWcode) | 0 | COWcode, 
           data = df)
m2 <- felm(protest ~ shutdown*event_time + I(case_lag/1000) + death_rate_lag + protest_lag | factor(COWcode) | 0 | COWcode, 
           data = df)
m3 <- felm(protest ~ shutdown*event_time + shutdown*I(event_time^2)  + I(case_lag/1000) + death_rate_lag + protest_lag| factor(COWcode) | 0 | COWcode, 
           data = df)
m4 <- felm(protest_bin ~ shutdown*event_time | factor(COWcode) | 0 | COWcode, 
           data = df)
m5 <- felm(protest_bin ~ shutdown*event_time + I(case_lag/1000) + death_rate_lag + protest_lag| factor(COWcode) | 0 | COWcode, 
           data = df)
m6 <- felm(protest_bin ~ shutdown*event_time + shutdown*I(event_time^2)  + I(case_lag/1000) + death_rate_lag + protest_lag | factor(COWcode) | 0 | COWcode, 
           data = df)

model1Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m1$beta[1],
                          SE = m1$cse[1],
                          Models = "Baseline")
model2Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m2$beta[1],
                          SE = m2$cse[1],
                          Models = "Covariates")
model3Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m3$beta[1],
                          SE = m3$cse[1],
                          Models = "Quadradic")
allModelFrame <- data.frame(rbind(model1Frame, model2Frame, model3Frame))  # etc.

interval1 <- -qnorm((1-0.9)/2)  
interval2 <- -qnorm((1-0.95)/2)  

# Plot
zp1 <- ggplot(allModelFrame, aes(shape = Models))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
zp1 <- zp1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2), fill = "WHITE")
zp1 <- zp1 + theme_tufte()
zp1 <- zp1 + ggtitle("Protest (Levels)")  + xlab("") + ylab("Shutdown Coefficient") + 
  scale_fill_discrete(name = 'Model') + theme(legend.position = "bottom")
zp1

model1Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m4$beta[1],
                          SE = m4$cse[1],
                          Models = "Baseline")
model2Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m5$beta[1],
                          SE = m5$cse[1],
                          Models = "Covariates")
model3Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m6$beta[1],
                          SE = m6$cse[1],
                          Models = "Quadradic")
allModelFrame <- data.frame(rbind(model1Frame, model2Frame, model3Frame))  # etc.

interval1 <- -qnorm((1-0.9)/2) 
interval2 <- -qnorm((1-0.95)/2)  

zp2 <- ggplot(allModelFrame, aes(shape = Models))
zp2 <- zp2 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
zp2 <- zp2 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp2 <- zp2 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2), fill = "WHITE")
zp2 <- zp2 + theme_tufte()
zp2 <- zp2 + ggtitle("Protests (Binary)") + xlab("") + ylab("Shutdown Coefficient") + 
  scale_fill_discrete(name = 'Model') + theme(legend.position = "bottom")

m1 <- felm(riot ~ shutdown*event_time | factor(COWcode) | 0 | COWcode, 
           data = df)
m2 <- felm(riot ~ shutdown*event_time + I(case_lag/1000) + death_rate_lag + riot_lag | factor(COWcode) | 0 | COWcode, 
           data = df)
m3 <- felm(riot ~ shutdown*event_time + shutdown*I(event_time^2)  + I(case_lag/1000) + death_rate_lag + riot_lag| factor(COWcode) | 0 | COWcode, 
           data = df)
m4 <- felm(riot_bin ~ shutdown*event_time | factor(COWcode) | 0 | COWcode, 
           data = df)
m5 <- felm(riot_bin ~ shutdown*event_time + I(case_lag/1000) + death_rate_lag + riot_lag | factor(COWcode) | 0 | COWcode, 
           data = df)
m6 <- felm(riot_bin ~ shutdown*event_time + shutdown*I(event_time^2)  + I(case_lag/1000) + death_rate_lag + riot_lag | factor(COWcode) | 0 | COWcode, 
           data = df)

model1Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m1$beta[1],
                          SE = m1$cse[1],
                          Models = "Baseline")
model2Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m2$beta[1],
                          SE = m2$cse[1],
                          Models = "Covariates")
model3Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m3$beta[1],
                          SE = m3$cse[1],
                          Models = "Quadradic")
allModelFrame <- data.frame(rbind(model1Frame, model2Frame, model3Frame))  

interval1 <- -qnorm((1-0.9)/2)  
interval2 <- -qnorm((1-0.95)/2)  

zp3 <- ggplot(allModelFrame, aes(shape = Models))
zp3 <- zp3 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
zp3 <- zp3 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp3 <- zp3 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2), fill = "WHITE")
zp3 <- zp3 + theme_tufte()
zp3 <- zp3 + ggtitle("Riots (Levels)")  + xlab("") + ylab("Shutdown Coefficient") + 
  scale_fill_discrete(name = 'Model') + theme(legend.position = "bottom")

model1Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m4$beta[1],
                          SE = m4$cse[1],
                          Models = "Baseline")
model2Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m5$beta[1],
                          SE = m5$cse[1],
                          Models = "Covariates")
model3Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m6$beta[1],
                          SE = m6$cse[1],
                          Models = "Quadradic")
allModelFrame <- data.frame(rbind(model1Frame, model2Frame, model3Frame))  # etc.

interval1 <- -qnorm((1-0.9)/2)  
interval2 <- -qnorm((1-0.95)/2) 

zp4 <- ggplot(allModelFrame, aes(shape = Models))
zp4 <- zp4 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
zp4 <- zp4 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp4 <- zp4 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2), fill = "WHITE")
zp4 <- zp4 + theme_tufte()
zp4 <- zp4 + ggtitle("Riots (Binary)") + xlab("") + ylab("Shutdown Coefficient") + 
  scale_fill_discrete(name = 'Model') + theme(legend.position = "bottom")

plot_2 <- cowplot::plot_grid(zp1, zp2, zp3, zp4, labels = 'AUTO')
ggsave('protest_plot.png', plot_2, width = 13)


#### Appendix Table of Countries #### 
df_check <- df %>% 
  group_by(cname) %>% 
  mutate(repress_pre = mean(repress[shutdown==0]), 
         repress_post = mean(repress[shutdown==1]))
df_check$change_repress <- round(df_check$repress_post-df_check$repress_pre, 1)

df_check <- df_check %>% 
  group_by(cname) %>% 
  summarise(repress_pre = round(mean(repress_pre),1), 
            repress_post = round(mean(repress_post), 1), 
            change_repress = round(mean(change_repress),1))

xtable::xtable(df_check)


#### Stringency Index (Continuous Treatment) #### 


m1 <- felm(repress ~ StringencyIndex*event_time | factor(COWcode) | 0 | COWcode, 
           data = df)
m2 <- felm(repress ~ StringencyIndex*event_time + I(case_lag/1000) + death_rate_lag + repress_lag | factor(COWcode) | 0 | COWcode, 
           data = df)
m3 <- felm(repress ~ StringencyIndex*event_time + StringencyIndex*I(event_time^2)  + I(case_lag/1000) + death_rate_lag + repress_lag| factor(COWcode) | 0 | COWcode, 
           data = df)
m4 <- felm(repress_bin ~ StringencyIndex*event_time | factor(COWcode) | 0 | COWcode, 
           data = df)
m5 <- felm(repress_bin ~ StringencyIndex*event_time + I(case_lag/1000) + death_rate_lag + repress_lag | factor(COWcode) | 0 | COWcode, 
           data = df)
m6 <- felm(repress_bin ~ StringencyIndex*event_time + StringencyIndex*I(event_time^2)  + I(case_lag/1000) + death_rate_lag + repress_lag | factor(COWcode) | 0 | COWcode, 
           data = df)

model1Frame <- data.frame(Variable = 'StringencyIndex',
                          Coefficient = m1$beta[1],
                          SE = m1$cse[1],
                          Model = "Baseline")
model2Frame <- data.frame(Variable = 'StringencyIndex',
                          Coefficient = m2$beta[1],
                          SE = m2$cse[1],
                          Model = "Covariates")
model3Frame <- data.frame(Variable = 'StringencyIndex',
                          Coefficient = m3$beta[1],
                          SE = m3$cse[1],
                          Model = "Quadradic")
allModelFrame <- data.frame(rbind(model1Frame, model2Frame, model3Frame)) 

interval1 <- -qnorm((1-0.9)/2)  
interval2 <- -qnorm((1-0.95)/2)  

zp1 <- ggplot(allModelFrame, aes(shape = Model))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
zp1 <- zp1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2), fill = "WHITE")
zp1 <- zp1 + theme_tufte()
zp1 <- zp1 + ggtitle("Repression (Levels)")  + xlab("") + 
  scale_fill_discrete(name = 'Model') + theme(legend.position = "bottom")
zp1


model1Frame <- data.frame(Variable = 'StringencyIndex',
                          Coefficient = m4$beta[1],
                          SE = m4$cse[1],
                          Model = "Baseline")
model2Frame <- data.frame(Variable = 'StringencyIndex',
                          Coefficient = m5$beta[1],
                          SE = m5$cse[1],
                          Model = "Covariates")
model3Frame <- data.frame(Variable = 'StringencyIndex',
                          Coefficient = m6$beta[1],
                          SE = m6$cse[1],
                          Model = "Quadradic")
allModelFrame <- data.frame(rbind(model1Frame, model2Frame, model3Frame))  

interval1 <- -qnorm((1-0.9)/2)  
interval2 <- -qnorm((1-0.95)/2)  

zp2 <- ggplot(allModelFrame, aes(shape = Model))
zp2 <- zp2 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
zp2 <- zp2 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp2 <- zp2 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2), fill = "WHITE")
zp2 <- zp2 + theme_tufte()
zp2 <- zp2 + ggtitle("Repression (Binary)") + xlab("") + 
  scale_fill_discrete(name = 'Model') + theme(legend.position = "bottom")

plot_1 <- plot_grid(zp1, zp2, labels = 'AUTO')

ggsave('string_plot.png', plot_1, width = 13)



#### Seasonal Adjustment #### 
df_sub <- subset(df, as.numeric(month) < 7 & as.numeric(month)>2)

m1 <- felm(repress ~ shutdown*event_time | factor(COWcode) | 0 | COWcode, 
           data = df_sub)
m2 <- felm(repress ~ shutdown*event_time + I(case_lag/1000) + death_rate_lag + repress_lag | factor(COWcode) | 0 | COWcode, 
           data = df_sub)
m3 <- felm(repress ~ shutdown*event_time + shutdown*I(event_time^2)  + I(case_lag/1000) + death_rate_lag + repress_lag| factor(COWcode) | 0 | COWcode, 
           data = df_sub)
m4 <- felm(repress_bin ~ shutdown*event_time | factor(COWcode) | 0 | COWcode, 
           data = df_sub)
m5 <- felm(repress_bin ~ shutdown*event_time + I(case_lag/1000) + death_rate_lag + repress_lag | factor(COWcode) | 0 | COWcode, 
           data = df_sub)
m6 <- felm(repress_bin ~ shutdown*event_time + shutdown*I(event_time^2)  + I(case_lag/1000) + death_rate_lag + repress_lag | factor(COWcode) | 0 | COWcode, 
           data = df_sub)
model1Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m1$beta[1],
                          SE = m1$cse[1],
                          Models = "Baseline")
model2Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m2$beta[1],
                          SE = m2$cse[1],
                          Models = "Covariates")
model3Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m3$beta[1],
                          SE = m3$cse[1],
                          Models = "Quadradic")

allModelFrame <- data.frame(rbind(model1Frame, model2Frame, model3Frame))  

interval1 <- -qnorm((1-0.9)/2) 
interval2 <- -qnorm((1-0.95)/2)  

zp1 <- ggplot(allModelFrame, aes(shape = Models))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
zp1 <- zp1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2), fill = "WHITE")
zp1 <- zp1 + theme_tufte()
zp1 <- zp1 + ggtitle("Repression (Levels)")  + xlab("") + ylab("Shutdown Coefficient") + 
  scale_fill_discrete(name = 'Model') + theme(legend.position = "bottom")
zp1


model1Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m4$beta[1],
                          SE = m4$cse[1],
                          Models = "Baseline")
model2Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m5$beta[1],
                          SE = m5$cse[1],
                          Models = "Covariates")
model3Frame <- data.frame(Variable = 'Shutdown',
                          Coefficient = m6$beta[1],
                          SE = m6$cse[1],
                          Models = "Quadradic")
allModelFrame <- data.frame(rbind(model1Frame, model2Frame, model3Frame))  

interval1 <- -qnorm((1-0.9)/2)  
interval2 <- -qnorm((1-0.95)/2)  

zp2 <- ggplot(allModelFrame, aes(shape = Models))
zp2 <- zp2 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
zp2 <- zp2 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp2 <- zp2 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2), fill = "WHITE")
zp2 <- zp2 + theme_tufte()
zp2 <- zp2 + ggtitle("Repression (Binary)") + xlab("") + ylab("Shutdown Coefficient") + 
  scale_fill_discrete(name = 'Model') + theme(legend.position = "bottom")

plot_1 <- cowplot::plot_grid(zp1, zp2, labels = 'AUTO')

ggsave('season_plot.png', plot_1, width = 13)




